Namibian hake model update, 2024

Author

John Kathena, Jim Ianelli

Published

July 16, 2024

Show the code
library(knitr)
hook_output = knit_hooks$get('output')
knit_hooks$set(output = function(x, options) {
  # this hook is used only when the linewidth option is not NULL
  if (!is.null(n <- options$linewidth)) {
    x = knitr:::split_lines(x)
    # any lines wider than n should be wrapped
    if (any(nchar(x) > n)) x = strwrap(x, width = n)
    x = paste(x, collapse = '\n')
  }
  hook_output(x, options)
})

knitr::opts_chunk$set(collapse = TRUE, comment = "  " ,fig.align = 'center', cache=FALSE,tidy.opts=list(width.cutoff=80), tidy=TRUE )
knitr::opts_knit$set(root.dir=here::here())
#knitr::opts_chunk$set(warning=F, message=F, echo=F, results=F,fig.width=6, fig.height=5)

Project overview

The following reflects my interpretation of the Marine Stewardship Council’s certification request. There is a need to catch up on missed milestones and outlines the necessary steps for the upcoming Year 4 milestone.

Key Points

  1. Year 3 Milestone Missed: The milestone required revised stock assessments for M. paradoxus, which was not met.
  2. Year 4 Milestone: By February 2025, the MFMR must use the Harvest Control Rule (HCR) systematically to verify the TAC for M. paradoxus. This needs to be applied in the August/September 2024 management meetings.
  3. Namibian Stock Assessment: There’s a recommendation to review and re-evaluate the assumptions and parameter values of assessment models, particularly the pessimistic base case model.
  4. Implementation Issues: MFMR has Dr. Ianelli’s report but not the code to run the model, and training is required for the Namibian team.

Draft agenda

The following draft agenda outlines the steps to address the missed milestones and prepare for the upcoming Year 4 milestone.

Week 1:

  1. Day 1-2: Review and Planning
    • Review the Year 3 milestone requirements and current progress.
    • Plan steps to implement the HCR for M. paradoxus.
  2. Day 3-4: Data Preparation
    • Gather and prepare Namibia stock assessment data.
    • Coordinate with MFMR to understand current data handling and management practices.
  3. Day 5: Meeting Preparation
    • Prepare documentation and a presentation for the MFMR management meeting.
    • Outline the steps needed for the August/September 2024 meeting to include HCR in TAC setting.

Week 2:

  1. Day 1-2: Model Review
    • Review Dr. Ianelli’s report and identify key elements for implementation.
    • Develop a preliminary implementation plan for the HCR model.
  2. Day 3-4: Training Coordination
    • Arrange a training session with Dr. Ianelli or another suitable individual for MFMR.
    • Coordinate with the training provider and MFMR to schedule the session.
  3. Day 5: Reporting
    • Compile a progress report summarizing activities, challenges, and next steps.
    • Send the report to Hugh and relevant stakeholders for feedback.

This agenda ensures a systematic approach to address the milestones and prepare for the upcoming management meeting, focusing on implementing the HCR and providing necessary training to MFMR.

Below are two main sections, first on model developments and second on application of a control rule that accounts for the signals in the data on the different species.

Assessment model runs

The original base-case model was evaluated for a number of features and extensions. These included focus on what data components were fit well and how improvements in consistency can be made. For the latter part, we found that the fits to the index and CPUE data were particularly poor and could be improved.

Below are examples which show how models can be run based on the directory location.

Model descriptions

The following table was developed based on testing the model with different assumptions and data sources. Key differences from the 2023 assessment configuration was the assumption that model estimation of variance terms was appropriate. This feature resulted in unacceptable residual patterns and essentially a complete down weighting of the index data. We used the assumed variance terms (CVs) for the indices in all of the following model configurations:

Model Description
Previous base case As specified in past assessments, estimated steepness and all variance terms
Base case (m0) Model with survey “minus group” to be ages 0, and 1 instead of 0, 1, and 2 as done in the past, steepness fixed at 0.7, q estimated, and time-varying fishery asymptotic selectivity specified.
m1 As base case but with survey catchability fixed at 1.0
m2 As base case but with survey catchability fixed at 0.5
m3 As base case but with natural mortality estimated
m4 As base case but with fishery selectivity allowed to be dome-shaped
m5 As base case but with stock-recruit steepness fixed at 0.43
m6 As base case but with stock-recruit steepness fixed at 0.9
Show the code
# Run names
mod_ref <- c("old_bc", "m0", "m1", "m2", "m3", "m4", "m5", "m6")
mod_dir <- c("old_bc", "m0", "m1", "m2", "m3", "m4", "m5", "m6")
mod_label <- c("2023 base case", "2024 base case", "Model 1", "Model 2", "Model 3",
    "Model 4", "Model 5", "Model 6")
res <- get_results(mod_names. = mod_label, moddir = mod_dir)
modlst <- res$modlst
old_bc <- modlst[[1]]
m0 <- modlst[[2]]
moddiag <- res$moddiag
dfsrr <- data.frame()
for (i in 1:length(mod_ref)) {
    # assign(modlst[[i]], run_nh(mod_dir[i] ,runit=FALSE) )
    dfsrr <- rbind(dfsrr, data.frame(Model = mod_label[i], SSB = modlst[[i]]$SSB,
        R = modlst[[i]]$Pred_Rec))
}
mods <- data.frame()
for (i in 1:length(mod_ref)) {
    mods <- rbind(mods, data.frame(moddiag[[i]], Model = names(moddiag[i])))
}
# head(mods)

Fits to index data

Base case model fits to main survey data compared to the previous assessment.

Base case model fits to main survey data compared to the previous assessment.

Base case model fits to the CPUE index 3 data compared to the previous assessment.

Base case model fits to the CPUE index 3 data compared to the previous assessment.
Show the code
dfcpue <- rbind(data.frame(Year = 1964:2023, Obs = old_bc$Obs_CPUE_1, predicted = old_bc$e_CPUE_1,
    Model = "Base case 2023"), data.frame(Year = 1964:2024, Obs = m0$Obs_CPUE_1,
    predicted = m0$e_CPUE_1, Model = "Base case 2024"))
dfcpue |>
    filter(Obs > 0) |>
    ggplot(aes(x = Year, y = Obs, color = Model)) + geom_point(color = "black") +
    geom_line(aes(y = predicted)) + geom_point(aes(y = predicted)) + ggtitle("CPUE 1") +
    ylim(0, NA)

Base case model fits to the CPUE index 1 data compared to the previous assessment.

Base case model fits to the CPUE index 1 data compared to the previous assessment.

SSB results

Base case model showing the SSB estimates compared to the previous assessment.

Base case model showing the SSB estimates compared to the previous assessment.
Show the code
mods |>
    filter(Model %in% mod_label[1:2], Year > 1960, Variable == "Depletion") |>
    ggplot(aes(x = Year, y = value, ymin = ymin, ymax = ymax, type = Model, fill = Model)) +
    geom_ribbon(alpha = 0.24) + ggthemes::theme_few() + geom_line(aes(color = Model)) +
    geom_point(aes(color = Model, shape = Model), size = 1) + ylab("Relative SSB") +
    xlab("Year")

Base case model showing the relative SSB estimates compared to the previous assessment.

Base case model showing the relative SSB estimates compared to the previous assessment.

Base case model showing the relative SSB estimates compared to the previous assessment.

Base case model showing the relative SSB estimates compared to the previous assessment.

Base case model showing the relative SSB estimates compared to the previous assessment.

Base case model showing the relative SSB estimates compared to the previous assessment.

Base case model showing the relative SSB estimates compared to the previous assessment.

Base case model showing the relative SSB estimates compared to the previous assessment.

Age composition fits

Base case model fits to survey and fishery age composition data. Note that the base case model uses a ‘minus group’ equal to ‘1’ for the survey data.

Base case model fits to survey and fishery age composition data. Note that the base case model uses a ‘minus group’ equal to ‘1’ for the survey data.

Base case model fits to survey and fishery age composition data. Note that the base case model uses a ‘minus group’ equal to ‘1’ for the survey data.

Base case model fits to survey and fishery age composition data. Note that the base case model uses a ‘minus group’ equal to ‘1’ for the survey data.

Selectivity

Selectivity estimates for the base-case model run.

Selectivity estimates for the base-case model run.

Stock-recruitment curves

The following plot shows the stock-recruitment curves for the base case and models 1, 2, 3, and 5 and 6.

   Warning: The shape palette can deal with a maximum of 6 discrete values because more
   than 6 becomes difficult to discriminate
   ℹ you have requested 8 values. Consider specifying shapes manually if you need
     that many have them.
   Warning: Removed 40 rows containing missing values or values outside the scale range
   (`geom_point()`).

Model runs

Model runs

Stock status comparisons

Unclear why this table doesn’t show up…

Show the code
dftmp <- NULL
mod_scen <- c(2:8)
filler <- " "
names(filler) <- "-ln(Likelihood)"
for (ii in mod_scen) {
    x <- modlst[[ii]]
    nll <- x$ObjFun
    names(nll) <- "Overall"
    CPUE <- x$CPUE_Like
    names(CPUE) <- "CPUE"
    surv <- x$Survey_Like
    names(surv) <- "Survey"
    caa <- x$CAA_Likelihood
    names(caa) <- "Commercial CAA"
    caas <- x$CAAS_Likelihood
    names(caas) <- "Survey CAA"
    rec <- x$RecRes_Likelihood
    names(rec) <- "Rec. resids."
    oneyrold <- x$Oneyearold_Likelihood
    names(oneyrold) <- "One yr-old biomass"
    NumPars <- x$Npars
    names(NumPars) <- "Number parameters"
    AIC <- x$Akaike
    names(AIC) <- "AIC"
    v <- c(filler, nll, CPUE, surv, caa, caas, rec, oneyrold, NumPars, AIC)
    dftmp <- cbind(dftmp, v)
}
dftmp <- data.frame(rownames(dftmp), dftmp, row.names = NULL)
names(dftmp) <- c("Component", mod_label[mod_scen])
tabcap <- "Model fits to data components"
tab <- xtable(dftmp, caption = tabcap, type = "html", label = "tab:mgt_quants", digits = 0)
# align=paste0('lll',strrep('r',length(mod_scen+1))))
kable(tab, caption.placement = "top", include.rownames = FALSE, sanitize.text.function = function(x) {
    x
})
Component 2024 base case Model 1 Model 2 Model 3 Model 4 Model 5 Model 6
-ln(Likelihood)
Overall -58.9499 -24.2444 31.237 -40.5486 -24.7886 -56.6279 -22.494
CPUE -46.289 -47.4902 -44.8216 -46.5351 -47.6547 -46.103 -46.5013
Survey -31.4817 -20.7397 -18.0356 -31.4318 -26.5236 -32.0501 -30.9079
Commercial CAA -125.707 -116.119 -72.0791 -123.738 -126.094 -123.984 -118.683
Survey CAA 93.2688 96.4544 119.897 91.1227 114.223 92.7225 92.6993
Rec. resids. 12.738 27.5852 18.2565 29.2642 10.9395 13.4603 41.5074
One yr-old biomass 68.1504 68.6746 71.4361 68.4665 69.1887 68.1901 68.4218
Number parameters 83 83 83 84 91 83 83
AIC 48.1002 117.511 228.474 86.9029 132.423 52.7442 121.012
Show the code
dftmp <- NULL
mod_scen <- c(2:8)
for (ii in mod_scen) {
    x <- modlst[[ii]]
    "-ln(Likelihood)"
    nll <- x$ObjFun
    names(nll) <- "Overall"
    CPUE <- x$CPUE_Like
    names(CPUE) <- "CPUE"
    surv <- x$Survey_Like
    names(surv) <- "Survey"
    caa <- x$CAA_Likelihood
    names(caa) <- "Commercial CAA"
    caas <- x$CAAS_Likelihood
    names(caas) <- "Survey CAA"
    rec <- x$RecRes_Likelihood
    names(rec) <- "Rec. resids."
    oneyrold <- x$Oneyearold_Likelihood
    names(oneyrold) <- "One yr-old biomass"
    NumPars <- x$Npars
    names(NumPars) <- "Number parameters"
    AIC <- x$Akaike
    names(AIC) <- "AIC"
    v <- c(nll, CPUE, surv, caa, caas, rec, oneyrold, NumPars, AIC)
    dftmp <- cbind(dftmp, v)
}
dftmp <- data.frame(rownames(dftmp), dftmp, row.names = NULL)
names(dftmp) <- c("Component", mod_label[mod_scen])
tabcap <- "Model fits to data components"
tab <- xtable(dftmp, caption = tabcap, type = "html", label = "tab:mgt_quants", digits = 0)
# align=paste0('lll',strrep('r',length(mod_scen+1))))
kable(tab, caption.placement = "top", include.rownames = FALSE, sanitize.text.function = function(x) {
    x
})
Component 2024 base case Model 1 Model 2 Model 3 Model 4 Model 5 Model 6
Overall -58.9499 -24.2444 31.2370 -40.5486 -24.7886 -56.6279 -22.4940
CPUE -46.2890 -47.4902 -44.8216 -46.5351 -47.6547 -46.1030 -46.5013
Survey -31.4817 -20.7397 -18.0356 -31.4318 -26.5236 -32.0501 -30.9079
Commercial CAA -125.7070 -116.1190 -72.0791 -123.7380 -126.0940 -123.9840 -118.6830
Survey CAA 93.2688 96.4544 119.8970 91.1227 114.2230 92.7225 92.6993
Rec. resids. 12.7380 27.5852 18.2565 29.2642 10.9395 13.4603 41.5074
One yr-old biomass 68.1504 68.6746 71.4361 68.4665 69.1887 68.1901 68.4218
Number parameters 83.0000 83.0000 83.0000 84.0000 91.0000 83.0000 83.0000
AIC 48.1002 117.5110 228.4740 86.9029 132.4230 52.7442 121.0120

But here the figure does… #| # x\(KspSTD # x\)KexpSTD
# Ksp # h # Bsp(2023) # Bexp(2023) # Bsp(2023)/Ksp # Bsp(2024) # Bexp(2024) # Bsp(2024)/Ksp # Bspmsy # Bspmsy/Ksp # Bsp(2023)/Bspmsy # Bsp(2024)/Bspmsy # MSY # AverageRY_5yrs # AverageRY_5yrs*Bsp/Bspmsy # ```{r mgt_quants, results = “asis”} # Mgt quants by model

Some reference point comparisons among model runs. aveRY_90 is the average replacement yeild since 1990, aveRY is the average replacement yield over the most recent 5 years before the current year, Cur_90 is the current (terminal year) SSB over the estimate from 1990, Cur_B0 is over the unfished estimate, and Cur_Bmsy is the ratio of current SSB over the Bmsy estimate.

Some reference point comparisons among model runs. aveRY_90 is the average replacement yeild since 1990, aveRY is the average replacement yield over the most recent 5 years before the current year, Cur_90 is the current (terminal year) SSB over the estimate from 1990, Cur_B0 is over the unfished estimate, and Cur_Bmsy is the ratio of current SSB over the Bmsy estimate.

Control rule application

Design of triggered control rule

The situation for developing a two-species control rule where catches between species and the trend in overall biomass for both species is combined is challenging. For management purposes, the goal is to avoid incidental takes in the proportion of one species that exceeds the historical levels of depletion for either stock. Fortunately, survey data are available that can be used to distinguish trends in the relative biomass for both stocks. The design of the triggered control rule therefore must consider patterns in the relative biomass from the survey data, the absolute biomass of the combined catch and biomass as modeled from the combined-stocks assessment. This provides a pragmatic approach using available data.

The steps in the control rule would be first to run a simple model that projected the survey biomass and relative proportion of M. paradoxus. Then, given the mean proportion over the period, compute the adjustment needed to the overarching control rule for the management procedure. For example, the historical range based on the survey has been between about one third of the mean value (in the earliest part of the period) to about 70% above the mean proportion. This range (especially the lower value) was used as a semi-empirical way to develop a minimum stock size threshold as part of the control rule. That is, when the stock of M. paradoxus drops below 30% of the mean proportion of the combined stocks estimate, the TAC recommendation for the combined stock would be zero. So if proportion of M. paradoxus \((p_y)\) is greater than 30% of the long-term mean proportion, then

\(TAC_y=(0.8 \sum_{y}^{y-4}\frac{RY_y} 5) \min(1,B_y/B_{MSY}) \min(1,p_y/\bar{p})\)

In words, the TAC in year y is equal to the catch under the current control rule times the ratio of the spawning biomass relative to BMSY (or proxy) and the externally estimated proportion of M. paradoxus from survey data. The second and third terms on the right hands side are depicted in Table 1 (and include the zero catch recommended if the M. paradoxus proportion dropped below the lowest value historically observed). The values of γ,λ were both set to 1.0 but could be adjusted (say to 0.5) to have a more gradual transition as the proportion or aggregate stock conditions drop below mean levels. We note that the specification of a survey linkage by the individual species provides a trigger (30% of long-term mean) that reduces the exploitation rate as the point of recruitment impairment (PRI) is approached.

Improvements to this approach could benefit from management that accounts for the depth of fishing Johnsen and Kathena (2012). These data were unavailable for inclusion in this study since they rely on predictions of the species mix based on the depth of fishing. Also, it is unclear that depth of fishing could be forecasted in a manner that could be built into a management measure.

The approached applied a smoother to each of the species time-series of survey data.

Modeling Namibian hake survey data by species

Fisheries stock assessments require data that are reliably collected and compiled. Secondarily, assessment models should be configured to match the assumptions associated with the observed data. To account for survey trends between the two species of hake we applied the estimated observation errors to a simple state-space random walk model. This approach has a number of options for how process-errors can be specified and estimated. The observation model applies the observation-error variances \((σ_(j,t)^2)\) for the \(j^{th}\) species in year \(t (x_(j,t) )\). The indices are fit to latent state variables, e.g., the underlying population trend \(ln(Z ̂_(j,t )\) as follows:

\(ln(Z_(j,t) ) = ln(Z ̂_(j,t) )+ϵ_(j,t),"where " ϵ_(j,t)∼N(0,σ_(j,t)^2 )\)

and the state equation and associated process error variance σ_PE^2 is defined as

\(ln(Z ̂_(j,t+1) )=ln(Z ̂_(j,t) )+η_(j,t),"where " η_(j,t)∼N(0,τ_j^2 ).\)

The process error variances τ_j^2 (which may or may not vary across indices) are fixed effect parameters and the unobserved species combined population ln(Z_(j,t) ) is estimated as a series of random effects. The model is fit using maximum likelihood estimation in TMB using the R package “rema” (Sullivan 2022). The survey data for each species was used with CVs applied for observation error specifications. The values for τ_j^2 were estimated for each species since their inherent variability may differ.

The above analysis provides a summary of the model runs and the design of a control rule that accounts for the signals in the data on the different species.

Notes

Throughout the weeks, we scrutinized data inputs and found a couple of inconsistencies. For example, data were provided for surveys in 2019 yet there were no observations in that year. Similarly, the abundance-at-length data from the surveys failed to identify strong periods of persistence by cohorts.

Ignore 7-vessel CPUE data set.